load packages and functions

source('mytheme.R')
# model version
modelversion = 'gap_log_rstan'
rstanmodelPath = 'modelrlt'
modelPath = paste0(rstanmodelPath, '/models/', modelversion)

Merge the model Result data

1 load all data and model results

AllExpData = read.csv(paste0("../data/AllValidData.csv")) %>%filter(valid == 1)
dur <- sort(unique(AllExpData$curDur))

AllExpData$WMSize <- factor(AllExpData$WMSize, labels = c("low", "medium",  "high")) 

# Replace first 3 chracters "Exp" with string "Exp. "
AllExpData$Exp <- gsub("^.{0,3}", "Exp. ", AllExpData$Exp)
AllExpData[which(AllExpData$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllExpData[which(AllExpData$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
## load model Prediction results
AllDat_predY <- read.csv(paste0(modelPath, "/rlt/AllDat_predY_",modelversion,".csv"))
AllDat_predY$WMSize <- as.factor(AllDat_predY$WMSize)
levels(AllDat_predY$WMSize) = c("low", "medium",  "high")
AllDat_predY$pred_Bias = AllDat_predY$mu_r - AllDat_predY$curDur
AllDat_predY$predErr = AllDat_predY$mu_r - AllDat_predY$repDur
AllDat_predY$relatErr = AllDat_predY$predErr / AllDat_predY$repDur

# rename experiments by replacing first 3 chracters "Exp" with string "Exp. "
AllDat_predY$Exp <- gsub("^.{0,3}", "Exp. ", AllDat_predY$Exp)
AllDat_predY[which(AllDat_predY$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllDat_predY[which(AllDat_predY$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
AllDat_predY$Exp = as.factor(AllDat_predY$Exp)

AllDat_predY$gap = 1
AllDat_predY[which(AllDat_predY$Gap == 2500),"gap"] = 2
AllDat_predY$gap <- factor(AllDat_predY$gap, labels = c("short", "long"))

2 Corrct rate

dat_Exp4b = dplyr::group_by(AllExpData%>%filter(Exp =='Exp. 4b'), WMSize, NSub, gap) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n = n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))%>%
  dplyr::group_by(gap, WMSize)%>%
  dplyr::summarize( n = n(),
                    mean_WMCrr = mean(m_WMCrr), se_WMCrr = sd(m_WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'WMSize', 'NSub'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'gap'. You can override using the `.groups` argument.
dat_Exp4b$gap<- factor(dat_Exp4b$gap, labels = c("short", "long")) 

plt_WMCrr_Exp4b = ggplot(data = dat_Exp4b, 
                         aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,group =gap, color = gap))+ 
  geom_line(aes(linetype = gap), stat = "identity",position = position_dodge(width = 0.3))+
  geom_point(aes(shape = gap), stat = "identity",position = position_dodge(width = 0.3))+ 
  geom_errorbar(width=.3,  position = position_dodge(width = 0.3)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  scale_color_manual(values=c("#999999", "#56B4E9"))+
  labs("Memory load", color = "Gap", shape = "Gap", linetype = "Gap", y = "Mean accuracy in WM task in Exp. 4b") +
  theme_new 
plt_WMCrr_Exp4b

model <-ezANOVA(data= AllExpData%>%filter(Exp =='Exp. 4b')%>%dplyr::group_by(NSub, WMSize, gap)%>%dplyr::summarise(m_WMCrr = mean(WMCrr)), dv = m_WMCrr, wid = NSub, within = .(WMSize, gap), type = 3, detailed = T)
## `summarise()` has grouped output by 'NSub', 'WMSize'. You can override using the
## `.groups` argument.
## Warning: Converting "NSub" to factor for ANOVA.
## Warning: "gap" will be treated as numeric.
## Warning: There is at least one numeric within variable, therefore aov() will be
## used for computation and no assumption checks will be obtained.
model
## $ANOVA
##       Effect DFn DFd          SSn        SSd          F            p p<.05
## 1     WMSize   2  30 5.670075e-01 0.10014543 84.9276173 4.426759e-13     *
## 2        gap   1  15 7.853375e-05 0.01091070  0.1079680 7.470119e-01      
## 3 WMSize:gap   2  30 1.273336e-03 0.03673923  0.5198813 5.998504e-01      
##           ges
## 1 0.793236200
## 2 0.000531086
## 3 0.008541940
model$ANOVA[4] / (model$ANOVA[4] + model$ANOVA[5])  
sub_WMCrr_Exp4b = AllExpData%>%filter(Exp =='Exp. 4b')%>%dplyr::group_by(NSub, WMSize, gap)%>%dplyr::summarise(mWMCrr = mean(WMCrr))%>%pivot_wider(names_from = c("WMSize", "gap"), values_from = c(mWMCrr), names_sep="_")
## `summarise()` has grouped output by 'NSub', 'WMSize'. You can override using the
## `.groups` argument.
write.csv(sub_WMCrr_Exp4b,paste0(modelPath, '/rlt/sub_WMCrr_Exp4b.csv'))
Exp.labs.1line <- c("Exp. 1 Control", "Exp. 2 Encoding", "Exp. 3 Reproduction", "Exp. 4a Both phases", "Exp. 4b Both phases\n \ \ \ \ \ \ \ \ \ \ \ \ \ with a gap")

plt_WMCrr <- ggplot(meanForPlot, aes(WMSize, mean_WMCrr, ymin = mean_WMCrr - se_WMCrr, ymax = mean_WMCrr + se_WMCrr,
                                     group =Exp, color = Exp))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  coord_cartesian(ylim = c(0.5, 1)) +
  labs(x = "Memory load", y = "Mean accuracy in WM task") +
  theme_new +scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9", "#1a9641", "#d7191c"), labels = Exp.labs.1line)+
  scale_shape_manual(values = Exp.labs.1line) 



plt_WMCrr_2<-ggarrange(plt_WMCrr, plt_WMCrr_Exp4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(5.5,4), labels = c("a", "b"))
plt_WMCrr_2

ggsave(paste0(getwd(), "/figures/plt_WMCrr_2.png"), plt_WMCrr_2, width = 9.5, height = 4)
### generate WM correct rates
AllExpData$WMCrr <- AllExpData$TPresent == AllExpData$WMRP
m_wmp<- dplyr::group_by(AllExpData, Exp, WMSize, NSub) %>% 
  dplyr::summarize(m_WMCrr = mean(WMCrr), n =n(), se_WMCrr = sd(WMCrr)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.

2.1 check model parameters

### Average Parameters 
mm_Baypar <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize( m_sig_s2 = mean(sig_s2),
                    m_sig_pr2_log = mean(sig_pr2_log),
                    m_ks= mean(ks), 
                    m_kr = mean(kr), 
                    m_ls = mean(ls), 
                    m_ts = mean(ts),
                    m_mu_pr = mean(mu_pr),
                    m_sig_pr2 = mean(sig_pr2),
                    m_mu_pr_log= mean(mu_pr_log),
                    m_sig_mn2 = mean(sig_mn2),
                    n= n(),
                    se_sig_s2 = sd(sig_s2)/sqrt(n-1),
                    se_sig_mn2 = sd(sig_mn2)/sqrt(n-1),
                    se_sig_pr2_log = sd(sig_pr2_log)/sqrt(n-1),
                    se_ks = sd(ks)/sqrt(n-1),
                    se_kr = sd(kr)/sqrt(n-1),
                    se_ls =sd(ls)/sqrt(n-1),
                    se_mu_pr_log = sd(mu_pr_log)/sqrt(n-1)) 
mm_Baypar

3 Prediction results

#calculate the mean reproduction biases for the five given intervals for all subjects
mpredY_sub <- dplyr::group_by(AllDat_predY, curDur, Exp, NSub, WMSize) %>%
  dplyr::summarize(n = n(), 
                   m_repDur = mean(repDur), 
                   sd_repDur = sd(repDur),
                   m_mu_r = mean(mu_r),
                   m_sig_r = mean(sig_r),
                   m_wp = mean(wp),
                   se_wp = sd(wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   cv =sd_repDur/ m_repDur,
                   pred_cv = mean(sig_r/mu_r),
                   predRP_err = mean(m_mu_r-m_repDur),
                   predVar_err = mean(m_sig_r-sd_repDur),
                   predRP_rerr = mean(abs(m_mu_r-m_repDur)/m_repDur),
                   predVar_rerr = mean(abs(m_sig_r-sd_repDur)/sd_repDur),
                   predcv_err = pred_cv-cv,
                   predcv_rerr = mean(abs(pred_cv-cv)/cv))
## `summarise()` has grouped output by 'curDur', 'Exp', 'NSub'. You can override
## using the `.groups` argument.
write_csv(dplyr::group_by(mpredY_sub, curDur, NSub) %>%
            dplyr::summarize(m_cv = mean(cv))%>%spread(curDur, m_cv), paste0(modelPath, '/rlt/m_cv.csv'))
## `summarise()` has grouped output by 'curDur'. You can override using the
## `.groups` argument.
mpredY_sub$RP_bias = mpredY_sub$m_repDur -mpredY_sub$curDur
mpredY_sub_new <- dplyr::group_by(mpredY_sub, curDur, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(curDur, m_RP_bias)
## `summarise()` has grouped output by 'curDur', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 1'), paste0(modelPath, '/rlt/RP_Bias_exp1.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 2'), paste0(modelPath, '/rlt/RP_Bias_exp2.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 3'), paste0(modelPath, '/rlt/RP_Bias_exp3.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 4a'), paste0(modelPath, '/rlt/RP_Bias_exp4a.csv'))
write_csv(mpredY_sub_new%>%filter(Exp == 'Exp. 4b'), paste0(modelPath, '/rlt/RP_Bias_exp4b.csv'))


mpredY_sub_WMsize <- dplyr::group_by(mpredY_sub, WMSize, Exp, NSub) %>%
  dplyr::summarize(m_RP_bias = mean(RP_bias))%>% spread(WMSize, m_RP_bias)
## `summarise()` has grouped output by 'WMSize', 'Exp'. You can override using the
## `.groups` argument.
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 3'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp3.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 4a'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4a.csv'))
write_csv(mpredY_sub_WMsize%>%filter(Exp == 'Exp. 4b'), paste0(modelPath, '/rlt/RP_Bias_WMsize_exp4b.csv'))
#### predicted data
m_predY <- mpredY_sub%>% 
  dplyr::group_by(Exp, curDur, WMSize) %>% 
  dplyr::summarize(m_m_repDur = mean(m_repDur),
                   m_sd_repDur = mean(sd_repDur),
                   m_m_sig_r =mean(m_sig_r),
                   m_m_mu_r = mean(m_mu_r),
                   m_m_wp = mean(m_wp),
                   n = n(),
                   m_se_wp = sd(se_wp)/sqrt(n-1),
                   log_lik =mean(log_lik),
                   mpredRP_err = mean(predRP_err),
                   mpredVar_err = mean(predVar_err),
                   mpredRP_rerr = mean(predRP_rerr),
                   mpredVar_rerr = mean(predVar_rerr),
                   cv= mean(cv),
                   pred_cv = mean(pred_cv),
                   mpredcv_err = mean(predcv_err),
                   mpredcv_rerr = mean(predcv_rerr))
m_predY_acc =  mpredY_sub%>% 
  dplyr::group_by(Exp) %>% 
  dplyr::summarize(mpred_rerr = mean(predRP_rerr)*100,
                   mpredVar_rerr = mean(predVar_rerr)*100, 
                   mpredcv_rerr = mean(predcv_rerr)*100)
m_predY_acc

4 WAIC and LOO-CV

#extract waic and loo-cv from parameter list
m_WAIC <- dplyr::group_by(Bayparlist, Exp) %>%
  dplyr::summarize(n =n(),
                   m_looic = mean(looic),
                   m_waic = mean(waic),
                   se_waic = sd(waic)/sqrt(n-1),
                   se_looic = sd(looic)/sqrt(n-1),
                   m_p_loo = mean(p_loo),
                   m_elpd_loo = mean(elpd_loo),
                   m_se_looic = mean(se_looic),
                   m_se_p_loo = mean(se_p_loo),
                   m_p_waic = mean(p_waic),
                   m_se_waic = mean(se_waic)) 

m_WAIC
#load test results
AllDat_newY <- read.csv(paste0(modelPath, "/rlt/AllDat_newY_",modelversion,".csv"))
AllDat_newY$WMSize <- as.factor(AllDat_newY$WMSize)
levels(AllDat_newY$WMSize) = c("low", "medium",  "high")
# Replace first 3 chracters "Exp" with string "Exp. "
AllDat_newY$Exp <- gsub("^.{0,3}", "Exp. ", AllDat_newY$Exp)
AllDat_newY[which(AllDat_newY$Exp == 'Exp. 4'),"Exp"] = "Exp. 4a"
AllDat_newY[which(AllDat_newY$Exp == 'Exp. 5'),"Exp"] = "Exp. 4b"
AllDat_newY$Exp = as.factor(AllDat_newY$Exp)

AllDat_newY$gap = 1
AllDat_newY[which(AllDat_newY$Gap == 2500, AllDat_newY$Exp == "Exp. 4b"), "gap"] = 2
AllDat_newY$gap = factor(AllDat_newY$gap, labels = c("short", "long"))

4.1 RP biase

Exp.labs.2lines <- c("Exp. 1\n Control", "Exp. 2\n Encoding", "Exp. 3\n Reproduction", "Exp. 4a\n Both", "Exp. 4b\n Both, with a gap")
names(Exp.labs.2lines) <- c("Exp. 1", "Exp. 2", "Exp. 3", "Exp. 4a", "Exp. 4b")
RP_bias  <- ggplot(data = m_predY%>%filter(Exp !="Exp. 4b"), aes(x = curDur, y = m_m_repDur - curDur, color=as.factor(WMSize), shape = as.factor(WMSize))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= m_newY%>%filter(Exp !="Exp. 4b"), aes(x=curDur, y=m_mu_r-curDur, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x=" ", y="Reproduction bias (s)", shape ="Memory Load", color = "Memory Load")+theme_new+
  colorSet3+
  theme(legend.position = "top")+ylim(-0.5, 0.4)

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias.png"), RP_bias, width = 6, height = 3)

RP_bias

RP_bias_Exp4b  <- ggplot(data = AllDat_predY %>% filter(Exp == 'Exp. 4b')%>% 
                           dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
                           dplyr::summarize(m_mu_r = mean(mu_r), m_repDur = mean(repDur), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x = curDur, y = m_repDur - curDur, group = interaction(gap, WMSize), color=as.factor(WMSize), shape = gap)) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= AllDat_newY %>% filter(Exp == 'Exp. 4b') %>% 
              dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
              dplyr::summarize(m_mu_r = mean(mu_r),  n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x=curDur, y=m_mu_r-curDur, group = interaction(WMSize, gap), linetype = gap, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x=" ", y="Reproduction bias (s)", shape ="Gap", linetype = "Gap")+theme_new+ guides(color = "none")+
  colorSet3+
  theme(legend.position = "top")+ylim(-0.5, 0.4)
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
RP_bias_Exp4b

## Figures in the MS
RP_bias_all<-ggarrange(RP_bias, RP_bias_Exp4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(4,1.4),  labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_all.png"), RP_bias_all, width = 9, height = 4.5)
RP_bias_all

RP_bias_1_4b <- ggplot(data = AllDat_predY %>% 
                         dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
                         dplyr::summarize(m_mu_r = mean(mu_r), m_repDur = mean(repDur), n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x = curDur, y = m_repDur - curDur, group = interaction(gap, WMSize), color=as.factor(WMSize), shape = as.factor(WMSize))) +
  geom_point(size=2, alpha = 0.5)+
  geom_line(data= AllDat_newY %>% 
              dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
              dplyr::summarize(m_mu_r = mean(mu_r),  n= n(), se_mu_r = sd(mu_r)/sqrt(n-1)), aes(x=curDur, y=m_mu_r-curDur, group = interaction(WMSize, gap), linetype = gap, color=WMSize)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x=" ", y="Reproduction bias (s)", shape ="Memory Load", linetype = "Gap", color = "Memory Load")+theme_new+
  colorSet3+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_1_4b.png"), RP_bias_1_4b, width = 9, height = 4)
RP_bias_1_4b

AllExpData[which(AllExpData$Exp != "Exp. 4b"), "gap"] = 1
AllExpData$gap = factor(AllExpData$gap, labels = c("short", "long"))

m_obsRP = AllExpData %>% dplyr::group_by(Exp, curDur, WMSize, NSub, gap) %>% 
                             dplyr::summarize(n = n(),
                                              m_repDur = mean(repDur),
                                              se_repDur = sd(repDur)/sqrt(n-1)) %>%
                             dplyr::group_by(Exp, curDur, WMSize, gap) %>% 
                             dplyr::summarize(m_m_repDur = mean(m_repDur),
                                              m_se_repDur = mean(se_repDur))
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize', 'NSub'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'curDur', 'WMSize'. You can override using the `.groups` argument.
RP_bias_obs_gap  <- ggplot(data = m_obsRP, aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize), linetype = factor(gap))) +
  geom_point()+
  geom_line()+
  geom_errorbar(width=.1,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur, linetype = gap)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape="Memory Load", color = "Memory Load", linetype = "Gap")+
  theme_new+colorSet3+
  scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")

ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs_gap.png"), RP_bias_obs_gap, width = 9, height = 4)

RP_bias_obs_gap

RP_bias_obs  <- ggplot(data = m_obsRP%>%filter(Exp != 'Exp. 4b'), aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize))) +
  geom_point()+
  geom_line()+
  geom_errorbar(width=.1,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", shape="Memory Load", color = "Memory Load")+
  theme_new+colorSet3+
  scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")

RP_bias_obs

RP_bias_obs_4b  <- ggplot(data = m_obsRP%>%filter(Exp == 'Exp. 4b'), aes(x = curDur, y = m_m_repDur-curDur,shape = as.factor(WMSize), color=as.factor(WMSize), linetype = factor(gap))) +
  geom_point()+
  geom_line()+
  geom_errorbar(width=.1,  aes(ymin = m_m_repDur-curDur - m_se_repDur, ymax = m_m_repDur -curDur + m_se_repDur, linetype = gap)) +
  geom_hline(yintercept = 0, linetype='dashed')+
  facet_grid(cols = vars(Exp), labeller = labeller(Exp = Exp.labs.2lines)) +
  labs(x="Sample intervals (s)", y="Reproduction bias(s)", linetype = "Gap")+
  theme_new+colorSet3+
  scale_x_continuous(breaks=seq(0, 1.6, 0.4))+ theme(legend.position="top")+ylim(-0.5, 0.2)+guides(color = "none") +guides(shape = "none")

RP_bias_obs_4b

## Figures in the MS
RP_bias_obs_all<-ggarrange(RP_bias_obs, RP_bias_obs_4b, common.legend = FALSE, ncol=2, nrow=1, widths = c(4,1.4),  labels = c("a", "b"))
ggsave(paste0(getwd(), "/", modelPath, "/figures/RP_bias_obs_all.png"), RP_bias_obs_all, width = 9, height = 4.5)
RP_bias_obs_all

5 weight of prior

plt_wp <- ggplot(data = AllDat_predY %>%dplyr::group_by(NSub, Exp, WMSize) %>% dplyr::summarise(m_wp = mean(wp)) %>%dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(mm_wp = mean(m_wp), n= n(), m_se_wp = sd(m_wp)/sqrt(n-1)), aes(Exp, mm_wp, ymin = mm_wp - m_se_wp, ymax = mm_wp + m_se_wp, group =interaction(Exp, WMSize), color = WMSize, shape = as.factor(WMSize)))+ 
  geom_line(stat = "identity",position = position_dodge(width = 0.2))+
  geom_point(stat = "identity",position = position_dodge(width = 0.2))+ 
  geom_errorbar(width=.2,  position = position_dodge(width = 0.2)) +
  colorSet5+
  labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load', shape = 'Memory Load') + #scale_x_discrete(labels= Exp.labs.2lines)+
  theme_new + theme(legend.position="top")
## `summarise()` has grouped output by 'NSub', 'Exp'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups` argument.
plt_wp
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_wp.png"), plt_wp, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
m_wp = AllDat_predY %>%dplyr::group_by(NSub, Exp, WMSize, gap) %>% dplyr::summarise(m_wp = mean(wp)) %>%dplyr::group_by(Exp, WMSize, gap) %>% dplyr::summarise(mm_wp = mean(m_wp), n= n(), m_se_wp = sd(m_wp)/sqrt(n-1))
## `summarise()` has grouped output by 'NSub', 'Exp', 'WMSize'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the `.groups` argument.
plt_wp_gap <- ggplot(data = m_wp, aes(Exp, mm_wp, ymin = mm_wp - m_se_wp, ymax = mm_wp + m_se_wp, group =interaction(Exp, WMSize, gap), color = WMSize, linetype = gap, shape = WMSize))+
  geom_line(stat = "identity", position = position_dodge(width = 0.3))+
  geom_point(stat = "identity",position = position_dodge(width = 0.3))+
  geom_errorbar(width=.3,  position = position_dodge(width = .3)) +
  colorSet5+
  labs(x = "", y = TeX("Weight of the prior $w_p$"), color = 'Memory Load', shape = 'Memory Load', linetype= 'Gap') +
  theme_new + theme(legend.position="top")

plt_wp_gap 
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

6 Indifference Point and slope

6.1 Inifference Point (Equation)

sub_wp = AllDat_predY %>%dplyr::group_by(NSub, Exp, gap) %>% dplyr::summarise(wp = mean(wp)) 
## `summarise()` has grouped output by 'NSub', 'Exp'. You can override using the
## `.groups` argument.
Bayparlist$gap = 'short'
Bayparlist_4b = Bayparlist%>%filter(Exp == 'Exp. 4b')
Bayparlist_4b$gap = 'long'
Bayparlist_all = rbind(Bayparlist, Bayparlist_4b)
# Indifference Point based on Equation 10 in manuscript
Eq_Inp_list <- dplyr::left_join(Bayparlist_all, sub_wp) %>% 
  dplyr::group_by(NSub, Exp)  %>%
  mutate(sig_post = sig_pr2_log * wp) %>%
  mutate(InP1_log = (kr*log(1)-(1-wp)*ks*log(1)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
  mutate(InP1 = exp(InP1_log)) %>%
  mutate(InP3_log = (kr*log(3)-(1-wp)*ks*log(3)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
  mutate(InP3 = exp(InP3_log)) %>%
  mutate(InP5_log = (kr*log(5)-(1-wp)*ks*log(5)+ sig_post/2 +wp*mu_pr_log)/wp)%>%
  mutate(InP5 = exp(InP5_log)) %>%
  dplyr::select(NSub, Exp,gap, InP1, InP3, InP5)
## Joining, by = c("NSub", "Exp", "gap")
write.csv(Eq_Inp_list, paste0(modelPath, '/rlt/Eq_Inp_list.csv'))
m_Eq_Inp_list = Eq_Inp_list %>%
  pivot_longer(cols = starts_with("InP"), names_to = "WMSize", names_prefix = "InP", values_to = "InP")
m_Eq_Inp_list$gap = factor(m_Eq_Inp_list$gap, levels = c("short", "long"))

m_Eq_Inp_list$WMSize <- factor(m_Eq_Inp_list$WMSize, labels = c("low", "medium",  "high")) 
plt_Inp_eq =  m_Eq_Inp_list%>%dplyr::group_by(Exp, WMSize, gap)%>%
  dplyr::summarise(m_inP= mean(InP), n =n(), se_inP = sd(InP)/sqrt(n()-1)) %>%ggplot(aes(x= Exp, y=m_inP, color = WMSize, shape = WMSize, linetype = factor(gap)))+
  geom_line(stat = "identity",position = position_dodge(width = 0.3))+
  geom_point(stat = "identity",position = position_dodge(width = 0.3))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.3)) +theme_new+
  labs(colour = "Memory Load", shape = "Memory Load", linetype = "Gap")+colorSet3+
  xlab(' ')+ylab("Indifference point (s)")+
  theme(legend.position = "top") 
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
plt_Inp_eq
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Inp_eq.svg"), plt_Inp_eq, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
dat = Eq_Inp_list%>%filter(Exp %in%c('Exp. 2', 'Exp. 3'))

dat$high_low_2 = 0.5 * (dat$InP5 - dat$InP1)
dat1 = dat %>% dplyr::select(c("Exp","NSub","high_low_2")) %>%pivot_wider(names_from =c('Exp'), values_from = c(high_low_2))

dat1$Exp2_3 = dat1$`Exp. 2` - dat1$`Exp. 3`
mean(dat1$Exp2_3)
## [1] -0.1490995
sd(dat1$Exp2_3)/sqrt(15)
## [1] 0.03349644
mean(dat1$`Exp. 3`)
## [1] 0.07411996
sd(dat1$`Exp. 3`)/sqrt(15)
## [1] 0.02502439

6.2 Inifference Point (curve)

AllDat_newY$predErr = AllDat_newY$mu_r - AllDat_newY$curDur

temp_newY <- AllDat_newY %>% filter(!(Exp == 'Exp. 4a' & Gap == 2500))  %>% select(Exp, WMSize, NSub, predErr, curDur, gap)

InP_curve<- temp_newY %>% dplyr::group_by(Exp, WMSize, NSub, gap)%>%
  dplyr::summarise(minErr = min(abs(predErr)), idx = which.min(abs(predErr)))
## `summarise()` has grouped output by 'Exp', 'WMSize', 'NSub'. You can override
## using the `.groups` argument.
InP_curve$InP_curve = temp_newY[InP_curve$idx,]$curDur 
InP_curve$y = temp_newY[InP_curve$idx,]$predErr + temp_newY[InP_curve$idx,]$curDur

#plot indifference points (the intersections of the Prediction curve with the diagonal)
plt_InP_curve<- ggplot(data = InP_curve%>%dplyr::group_by(Exp, WMSize, gap)%>% dplyr::summarise(m_InP = mean(InP_curve), se_InP = sd(InP_curve)/sqrt(n()-1)), aes(x= Exp, y=m_InP, color = as.factor(WMSize), linetype = factor(gap), shape = as.factor(WMSize)))+
  geom_line(stat = "identity",position = position_dodge(width = 0.3))+
  geom_point(stat = "identity",position = position_dodge(width = 0.3))+ 
  geom_errorbar(width=.3,  aes(ymin = m_InP - se_InP, ymax = m_InP + se_InP), position = position_dodge(width = 0.3)) +theme_new+
  labs(colour = "Memory Load", shape ="Memory Load")+colorSet3+
  xlab(' ')+ylab("indifference point (s)")+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_InP_curve.png"), plt_InP_curve, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_InP_curve
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

6.2.1 Observed RP

#Observed Indifference Point for Exp. 4b
obs_model <- function(df) {
  lm(repDur ~ curDur, data = df)
}
#Observed Indifference Point
obs_Inp_list <- AllDat_predY %>% 
  dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  dplyr::select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
## Warning: The `.drop` argument of `unnest()` is deprecated as of tidyr 1.0.0.
## All list-columns are now preserved.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
obs_Inp_list$model = NULL
obs_Inp_list$data = NULL
obs_Inp_list$inP = obs_Inp_list$Intercept /(1-obs_Inp_list$slope)
obs_Inp_list$slope = -1 * obs_Inp_list$slope
obs_Inp_list_no_gap <- AllDat_predY %>%
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%
  mutate(model = map(data, obs_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  dplyr::select(-std.error,-statistic, -p.value) %>%
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, slope = curDur)  # rename columns
obs_Inp_list_no_gap$model = NULL
obs_Inp_list_no_gap$data = NULL
obs_Inp_list_no_gap$inP = obs_Inp_list_no_gap$Intercept /(1-obs_Inp_list_no_gap$slope)
obs_Inp_list_no_gap$slope = -1 * obs_Inp_list_no_gap$slope
ezANOVA(data = obs_Inp_list%>%filter(Exp =='Exp. 4b'), dv= inP, wid=NSub, within= .(gap, WMSize) )
## Warning: Converting "NSub" to factor for ANOVA.
## $ANOVA
##       Effect DFn DFd        F           p p<.05         ges
## 2        gap   1  15 4.065468 0.062041074       0.019543220
## 3     WMSize   2  30 8.762649 0.001006806     * 0.078253848
## 4 gap:WMSize   2  30 0.584382 0.563670241       0.003061788
## 
## $`Mauchly's Test for Sphericity`
##       Effect         W           p p<.05
## 3     WMSize 0.4515805 0.003829536     *
## 4 gap:WMSize 0.5210587 0.010428128     *
## 
## $`Sphericity Corrections`
##       Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF] p[HF]<.05
## 3     WMSize 0.6458198 0.004991546         * 0.6808339 0.004255437         *
## 4 gap:WMSize 0.6761593 0.502620038           0.7194299 0.512209537
# plot the observed indifference points and slopes of RP
plt_obs_InP_slope_err<- ggplot(data = obs_Inp_list %>% group_by(Exp, WMSize)%>%
                                 dplyr::summarise(n=n(),
                                                  m_inP = mean(inP),
                                                  se_inP = sd(inP)/sqrt(n-1),
                                                  m_slope = mean(slope),
                                                  se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y=m_inP, color = WMSize, shape = WMSize))+
  geom_line(stat = "identity")+
  geom_point(stat = "identity")+ 
  geom_errorbar(width = 0.02, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP)) +
  geom_errorbarh(height =0.02, aes(xmin = m_slope - se_slope, xmax = m_slope + se_slope)) +
  theme_new+
  labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
  facet_grid(~Exp)+
  xlab('slope of reproduction bias')+ylab("indifference point (s)")+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_obs_InP_slope_err.png"), plt_obs_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_obs_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

6.2.2 Predicated RP

#Predicated Indifference Point for Exp.4b
pred_model <- function(df) {
  lm(mu_r ~ curDur, data = df)
}

pred_Inp_list <- AllDat_predY %>% 
  dplyr::group_by(NSub, Exp, WMSize, gap) %>% nest()  %>%
  mutate(model = map(data, pred_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  dplyr::select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, pred_slope = curDur)  # rename columns
pred_Inp_list$model = NULL
pred_Inp_list$data = NULL
pred_Inp_list$pred_inP = pred_Inp_list$Intercept /(1-pred_Inp_list$pred_slope)
pred_Inp_list$pred_slope = -1 * pred_Inp_list$pred_slope
pred_Inp_slope_no_gap <- AllDat_predY %>% 
  dplyr::group_by(NSub, Exp, WMSize) %>% nest()  %>%
  mutate(model = map(data, pred_model)) %>%  # linear regression
  mutate(slope = map(model, broom::tidy)) %>%  # get estimates
  unnest(slope, .drop = TRUE) %>% # remove raw data
  dplyr::select(-std.error,-statistic, -p.value) %>%  # remove unnessary columns
  spread(term, estimate) %>%   # spread stimates
  dplyr::rename(Intercept = `(Intercept)`, pred_slope = curDur)  # rename columns
pred_Inp_slope_no_gap$model = NULL
pred_Inp_slope_no_gap$data = NULL
pred_Inp_slope_no_gap$pred_inP = pred_Inp_slope_no_gap$Intercept /(1-pred_Inp_slope_no_gap$pred_slope)
pred_Inp_slope_no_gap$pred_slope = -1 * pred_Inp_slope_no_gap$pred_slope
m_pred_Inp_slope_no_gap = pred_Inp_slope_no_gap %>% group_by(Exp, WMSize)%>%
  dplyr::summarise(n=n(),
                   m_Intercept = mean(Intercept),
                   se_Intercept= sd(Intercept)/sqrt(n-1),
                   m_pred_inP = mean(pred_inP),
                   se_pred_inP = sd(pred_inP)/sqrt(n-1),
                   m_pred_slope = mean(pred_slope),
                   se_pred_slope = sd(pred_slope)/sqrt(n-1))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
# plot the observed indifference points and slopes of RP
plt_pred_InP_slope_err<- ggplot(data = m_pred_Inp_slope_no_gap, aes(x= m_pred_slope, y=m_pred_inP, color = WMSize, shape = WMSize))+
  geom_line(stat = "identity")+geom_point(stat = "identity")+ 
  geom_errorbar(width = 0.02, aes(ymin = m_pred_inP - se_pred_inP, ymax = m_pred_inP + se_pred_inP)) +
  geom_errorbarh(height =0.02, aes(xmin = m_pred_slope - se_pred_slope, xmax = m_pred_slope + se_pred_slope)) +
  geom_point(data = obs_Inp_list%>% group_by(Exp, WMSize)%>%
               dplyr::summarise(n=n(),
                                m_inP = mean(inP),
                                se_inP = sd(inP)/sqrt(n-1),
                                m_slope = mean(slope),
                                se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y =m_inP, color = WMSize, shape = WMSize))+
  theme_new+
  labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
  facet_grid(~Exp)+
  xlab('slope of reproduction')+ylab("indifference point (s)")+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_pred_InP_slope_err.png"), plt_pred_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_pred_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

InP_obs<-  ggplot(data = obs_Inp_list_no_gap %>%dplyr::group_by(WMSize, Exp) %>%dplyr::summarise(m_inP = mean(inP), se_inP = sd(inP)/sqrt(n()-1)), aes(x= Exp, y=m_inP, color = WMSize, shape = WMSize))+
  geom_line(stat = "identity",position = position_dodge(width = 0.3))+
  geom_point(stat = "identity",position = position_dodge(width = 0.3))+ 
  geom_errorbar(width=.3,  aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP), position = position_dodge(width = 0.3)) +theme_new+
  labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+
  xlab(' ')+ylab("observed indifference point (s)")+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'WMSize'. You can override using the
## `.groups` argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/InP_obs.png"), InP_obs, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_obs
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

InP_pred<-  ggplot(data = m_pred_Inp_slope_no_gap, aes(x= Exp, y=m_pred_inP, color = WMSize, shape = WMSize))+
  geom_line(stat = "identity",position = position_dodge(width = 0.3))+
  geom_point(stat = "identity",position = position_dodge(width = 0.3))+ 
  geom_errorbar(width=.3,  aes(ymin = m_pred_inP - se_pred_inP, ymax = m_pred_inP + se_pred_inP), position = position_dodge(width = 0.3)) +theme_new+
  labs(colour = "Memory Load", shape = "Memory Load")+colorSet3+ #scale_x_discrete(labels= Exp.labs.2lines)+
  xlab(' ')+ylab("indifference point (s)")+
  theme(legend.position = "top")


ggsave(paste0(getwd(), "/", modelPath, "/figures/InP_pred.png"), InP_pred, width = 3, height = 3)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
InP_pred
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

### Calculate predication error
Inp_list_no_gap = left_join(obs_Inp_list_no_gap, pred_Inp_slope_no_gap, by = c("NSub", "Exp", "WMSize"))
Inp_list_no_gap$InP_err = Inp_list_no_gap$pred_inP -Inp_list_no_gap$inP
Inp_list_no_gap$InP_rerr = 100*Inp_list_no_gap$InP_err/ Inp_list_no_gap$inP

Inp_list_no_gap$slope_err = Inp_list_no_gap$pred_slope - Inp_list_no_gap$slope
Inp_list_no_gap$slope_rerr = 100* Inp_list_no_gap$slope_err/Inp_list_no_gap$slope

m_Inp_list_no_gap = Inp_list_no_gap %>% dplyr::group_by(Exp) %>% dplyr::summarise(m_InP_rerr = mean(InP_rerr), m_slope_rerr = mean(slope_rerr), m_InP_rerr_abs = mean(abs(InP_rerr)), m_slope_rerr_abs = mean(abs(slope_rerr)))

m_Inp_list_no_gap$InP_auc = 100- m_Inp_list_no_gap$m_InP_rerr_abs
m_Inp_list_no_gap$slope_auc = 100-  m_Inp_list_no_gap$m_slope_rerr_abs

x

7 plot figures

#plot the predicated indifference points and slope of predicated RP
plt_pred_InP_slope_err<- ggplot(data = obs_Inp_list_no_gap%>% dplyr::group_by(Exp, WMSize)%>%
                                  dplyr::summarise(n=n(),
                                                   m_inP = mean(inP),
                                                   se_inP = sd(inP)/sqrt(n-1),
                                                   m_slope = mean(slope),
                                                   se_slope = sd(slope)/sqrt(n-1)), aes(x= m_slope, y=m_inP, shape = WMSize, color = WMSize))+
  geom_line(stat = "identity")+
  geom_point(stat = "identity")+ 
  geom_errorbar(width = 0.02, aes(ymin = m_inP - se_inP, ymax = m_inP + se_inP)) +
  geom_errorbarh(height =0.02, aes(xmin = m_slope - se_slope, xmax = m_slope + se_slope)) +
  theme_new+
  labs(colour = "Memory Load", shape = "Memory Load") +colorSet3+
  facet_grid(~Exp)+
  xlab('slope of reproduction')+ylab("indifference point (s)")+
  theme(legend.position = "top")
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_pred_InP_slope_err.png"), plt_pred_InP_slope_err, width = 5, height = 5)
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
plt_pred_InP_slope_err
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?

## Figures in the MS
fig3<-ggarrange(RP_bias, plt_pred_InP_slope_err, common.legend = TRUE, ncol=1, nrow=2,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig3.png"), fig3, width = 6, height = 5)
fig3

## combine InP and wp 
fig_wp_InP_gap<-ggarrange(plt_wp_gap, plt_Inp_eq, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig_wp_InP_gap.png"), fig_wp_InP_gap, width = 6, height = 3)
fig_wp_InP_gap

## combine InP and wp 
fig4<-ggarrange(plt_wp, InP_pred, common.legend = TRUE, ncol=2, nrow=1,  labels = c("a", "b"))
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
## geom_path: Each group consists of only one observation. Do you need to adjust
## the group aesthetic?
ggsave(paste0(getwd(), "/", modelPath, "/figures/fig4.png"), fig4, width = 6, height = 3)
fig4

8 Model prediction error

m_predErr_sub<- mpredY_sub%>% 
  dplyr::group_by(Exp, WMSize, NSub) %>% dplyr::summarise(
    mpredRP_err=mean(predRP_err),
    mpredVar_err=mean(predVar_err),
    mpredcv_err = mean(predcv_err),
    mpredRP_rerr = mean(predRP_rerr),
    mpredVar_rerr = mean(predVar_rerr),
    mpredcv_rerr = mean(predcv_rerr))
## `summarise()` has grouped output by 'Exp', 'WMSize'. You can override using the
## `.groups` argument.
m_predErr<- m_predY%>% 
  dplyr::group_by(Exp, WMSize) %>% dplyr::summarise(
    mmpredcv_err = mean(mpredcv_err),
    mmpredRP_err=mean(mpredRP_err),
    mmpredVar_err=mean(mpredVar_err),
    mmpredRP_rerr = mean(mpredRP_rerr),
    mmpredVar_rerr = mean(mpredVar_rerr),
    mmpredcv_rerr = mean(mpredcv_rerr))
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predErr

9 Model comparison (logarithmic vs. linear)

m_predErr_sub$model = 'logarithmic'
m_predErr$model = 'logarithmic'
linear_model = 'gap_linear_rstan'
m_predErr_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_", linear_model, ".csv"))
m_predErr_linear$X = NULL
m_predErr_sub_linear = read.csv(paste0(getwd(), "/", rstanmodelPath, '/models/', linear_model, "/rlt/m_predErr_sub_", linear_model, ".csv"))
m_predErr_sub_linear$X = NULL

m_predErr_sub_all = rbind(m_predErr_sub, m_predErr_sub_linear) 
m_predErr_all = rbind(m_predErr, m_predErr_linear)
m_predErr_all$WMSize = as.factor(m_predErr_all$WMSize)
levels(m_predErr_all$WMSize) =  c("low", "medium",  "high")
temp = m_predErr_all %>% filter(model == 'logarithmic') %>%summarise(abs_mmpredcv_err = abs(mmpredcv_err)) 
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
plt_Err_CV_all = ggplot(m_predErr_all, aes(abs(mmpredRP_err), abs(mmpredcv_err), group = interaction(model, Exp, WMSize), color = model, shape = Exp, size = WMSize)) + 
  geom_point() +
  geom_hline(yintercept = round(max(temp$abs_mmpredcv_err), 4)+0.0005, linetype='dashed')+
  xlab('Prediction error in the RP means (s)')+ ylab('Prediction error in CV')+colorSet3+
  scale_shape_manual(values = c(6, 7, 16, 17,8)) +
  theme_new+ 
  theme(legend.position = 'top')+
  labs(size = 'Memory Load')+
  guides(colour = guide_legend(order = 1, nrow=2,byrow=TRUE),
         shape = guide_legend(order =2, nrow=2,byrow=TRUE),
         size = guide_legend(order = 3, nrow=3,byrow=TRUE))

ggsave(paste0(getwd(), "/", modelPath, "/figures/plt_Err_CV_all.png"), plt_Err_CV_all, width = 7, height = 4)
## Warning: Using size for a discrete variable is not advised.
plt_Err_CV_all
## Warning: Using size for a discrete variable is not advised.

m_predY_acc =  m_predErr_sub_all%>% 
  dplyr::group_by(Exp, model) %>% 
  dplyr::summarize(mmpredRP_rerr = mean(mpredRP_rerr)*100,
                   mmpredVar_rerr = mean(mpredVar_rerr)*100,
                   mmpredcv_rerr = mean(mpredcv_rerr)*100,
                   mmpredRP_acc = (1-mean(mpredRP_rerr))*100,
                   mmpredVar_acc = (1-mean(mpredVar_rerr))*100,
                   mmpredCV_acc = (1-mean(mpredcv_rerr))*100)
## `summarise()` has grouped output by 'Exp'. You can override using the `.groups`
## argument.
m_predY_acc

Export data for spss